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We present a detailed analysis of the modulational insta- 
bility (MI) of ground-state bright solitary solutions of two in- 
coherently coupled nonlinear Schrodinger equations. Varying 
the relative strength of cross-phase and self-phase effects we 
show existence and origin of four branches of MI of the two- 
wave solitary solutions. We give a physical interpretation of 
our results in terms of the group velocity dispersion (GVD) 
induced polarization dynamics of spatial solitary waves. In 
particular, we show that in media with normal GVD spatial 
symmetry breaking changes to polarization symmetry break- 
ing when the relative strength of the cross-phase modulation 
exceeds a certain threshold value. The analytical and nu- 
merical stability analyses are fully supported by an extensive 
series of numerical simulations of the full model. 



I. INTRODUCTION 

The phenomenon of modulational instability (MI) can 
be defined as self-induced break up of an initially homo- 
geneous wave during its evolution in a nonlinear medium. 
Study of this phenomenon has been initiated in the 60 's, 
when MI was predicted in plasma physics M, nonlinear 
optics [g and physics of fluids [|| , and also observed ex- 
perimentally in form of filamentation of an optical beam 
propagating in an organic liquid B . Since that time MI 
has remained as one of the major topics of theoretical 
and experimental research in nonlinear physics and, in 
particular, in nonlinear physics of conservative systems 
|§|"§,0 . We will deal below with one classical example of 
such systems, unifying a number of previous results, and 
presenting new MI phenomena. Our approach stresses 
the central role of symmetries. 

A general formulation of the problem of nonlinear wave 
propagation via fundamental sets of equations, as for ex- 
ample, the Maxwell or Navier Stokes equations, is a very 
demanding task even for modern computers. Therefore a 
number of simplified models have been introduced which 
approximately describe either propagation of the wave 
itself, e.g., KdV equation jlj], or propagation of slowly 
varying wave envelope, e.g., the Nonlinear Schrodinger 
(NLS) equation |l|. 

The simplest solutions of envelope equations are con- 
tinuous wave (CW) solutions homogeneous in space and 
time. CWs in single NLS equation exhibit MI in cases 
when nonlinearity and group velocity dispersion (GVD) 



or diffraction act in opposition, e.g. when nonlinearity 
is positive GVD must be anomalous and if nonlinear- 
ity is negative GVD must be normal. This rule changes 
when, accounting for polarization, for different directions 
of wave vectors, or for different carrier frequencies of the 
interacting waves, one replaces the single NLS by the set 
of incoherently coupled NLS equations. Then, if non- 
linear coupling is strong enough, MI of CWs becomes 
possible for any signs of nonlinearity and GVD jl^-g^l ■ 

Another important class of solutions of nonlinear equa- 
tions are solitary solutions ('solitons'). They may also 
exhibit MI if they are localized in some dimensions but 
extended in one or more others. Mis of the envelope 
solitons of the single NLS and of the KdV solitons were 
pioneered, respectively, by Zakharov and Rubenchik |2l| ] 
and by Kadomtsev and Petviashvily p3|. Latter MI was 
studied in a number of other theoretical and experimen- 
tal works. For reviews on MI of bright and dark solitary 
waves see, respectively, p-[lO[ and |lO| , [ll| ]. 

From a formal point of view the problem of the solitary 
wave MI can be considered as continuation of the soliton 
spectrum at zero modulational frequency Q into the re- 
gion ri y^ 0. An important class of discrete eigenmodes 
at ri = are the zero eigenvalue (or neutral) modes. 
Presence of which is directly linked to symmetries of 
the model equations. On a qualitative level, similarities 
and differences between MI of solitons and CW solutions 
can be understood on the basis of a comparison between 
the corresponding neutral modes. For example, the ID, 
bright spatial soliton of NLS is modulationally unstable 
in media with either anomalous or normal GVD. In first 
case, the neutral mode associated with the phase sym- 
metry is excited ('neck' MI) and in latter situation the 
translational mode associated with a shift along the di- 
rection perpendicular to the wave propagation becomes 
unstable ('snake' MI) pit]. The phase mode is present as 
well for CW solution and this leads to MI for anomalous 
GVD. However, the translational mode of CW solution 
is null and therefore CWs are stable for normal GVD. 

Increasing the number of free parameters can lead to 
more complex scenarios of MI, because coexistence and 
competition between different types of instability are 
likely to happen. In a recent Letter |23| we considered 
GVD induced MI of spatial solitons due to nondegener- 
ate three- wave mixing. It was shown that presence of an 
additional phase symmetry leads to the appearance of a 
new branch of neck MI in media with normal GVD. It 
was found that the novel instability strongly dominates 



the usual snake MI throughout the region of sohton ex- 
istence. Note, that because of this dominance, physical 
mechanisms responsible for the relative strength of neck 
and snake instabilities remain to be understood. Among 
others, this last problem is addressed in the present work, 
where we study MI of the bright solitary solutions in the 
incoherently coupled NLS equations. The incoherent na- 
ture of the coupling results in the presence of two phase 
symmetries. In spite of the similar symmetry properties, 
the features of MI induced dynamics of solitons in the 
present model appear to be more rich compare to the 
three-wave mixing case. In particular we show that the 
relative strength of the nonlinear cross-coupling governs 
the competition between neck and snake Mis in media 
with normal GVD. 

The rest of this paper is organized as follows. In Sec- 
tion 2 we introduce model equations and discuss their 
physical relevance. In Section 3 the problem of MI of 
the solitary waves is formulated in general terms. MI of 
different kinds of solitary solutions and its physical inter- 
pretation in terms of polarization dynamics are detailed 
in Sections 4 and 5. Discussion and summary of main 
results are presented in Sections 6 and 7. 



II. APPLICATIONS OF INCOHERENTLY 

COUPLED NLS EQUATIONS TO PROPAGATION 

OF ELECTROMAGNETIC WAVES 

The evolution of two suitably scaled slowly varying in- 
coherently coupled wave envelopes Ei and E2 in a weakly 
nonlinear, dispersive and diffractive medium is governed 
by the following equations |13 



id.Ei + aiVlEi + jid^Ei 



{\E,\^ + P\E2\^)E,=Q, (1) 



id,E2 + a2V\E2 + j2dfE2 + (1^2 1' + P\Ei\'')E2 = 0, 

where Vj_ = idx + jdy. Longitudinal (z) and transverse 
(x, y) coordinates are respectively measured in units of 
a suitable diffraction length Idif and of a characteristic 
transverse size of the envelope. The coordinate t is the 
retarded time scaled to the parameter T^Jldif /Idis, where 
T is the temporal duration of the envelope and Idis is a 
characteristic GVD length. Diffraction parameters ai^2 
are positive while GVD parameters 71^2 can have either 
sign. Rescaling x,y,t once more one can always choose 
ai/a2 and I71I/I72I to be any convenient constants. The 
parameter f3 measures the relative strength of cross-phase 
modulation compare to self-phase modulation. The non- 
linearity was chosen to be self-focusing because below we 
are interested in the dynamics of bright solitary waves. 

Eqs. (1) describe a variety of physical situations but 
we will focus here on their application to propagation 
of electro-magnetic (e/m) waves. Using a circular polar- 
ization basis to describe propagation of quasimonochro- 
matic e/m waves in isotropic dielectric materials leads 
to Eqs. (1), where, in such a case, Ei and E2 are en- 
velopes of the left- and right- polarized components [E4| . 



The diffraction and GVD parameters can be taken as 

"1,2 = 0.5, 71,2 = 7 = ±0.5, and fH = 1 + Xxyyx/Xxxyy 

p3, where Xijii ^^ ^^^ nonlinear susceptibility tensor. 
For example, p = 2 for the nonrcsonant electronic non- 
linearity and (3 — 7 for the nonlinearity due to molec- 
ular orientation |24|. For e/m waves propagating in an 
isotropic plasma, /? strongly depends on the ratio between 
the frequency of the e/m wave and the characteristic elec- 
tron plasma frequency, and it can have either sign [[l5[ . 
Eqs. (1) can also be applied to describe interaction of cir- 
cularly polarised waves in waveguides filled with linearly 
isotropic material, as, e.g., CS2 liquid p5[ , in opposite 
situation linear coupling between Ei and E2 should be 
incorporated in the model. 

Counterpropagation of scalar waves in Kerr media 
obeys Eqs. (1) with f3 determined by the wavelength- 
scale refractive index gratings written by the interference 
pattern [0 . The value of /3 in this situation is directly 
linked with diffusion which washes out the grating mak- 
ing 1 < /? < 2 (/? = 2 for zero diffusion). Envelopes of 
incoherent copropagating waves in Kerr media also obey 
Eqs. (1) with (3 = 2 iQ. In these two situations the 
group velocity difference of the wave envelopes, which is 
not explicitly written in Eqs. (1), can be removed by a 
suitable phase shift. 

The limiting case f3 — > -t-00 describes a situation with 
zero self-modulation effects. This approximates the so 
called cascading limit of non-degenerate three- wave mix- 
ing in the quadratically nonlinear media ||26| . Therefore 
one can expect that for the large enough (3 MI of the 
solitary solutions of Eqs. (1) should be equivalent to the 
MI of the three-wave solitons p3] but that it should well 
be different for the relatively small (3. 

Because the discussion of a wide range of f3 values is 
more realistic in the context of the interaction of the cir- 
cularly polarized waves, below we mainly use terminology 
which is appropriate to this case. 



III. MODULATIONAL INSTABILITY OF 

SOLITONS. GENERAL FORMULATION OF THE 

PROBLEM 

The primary target of the present paper is understand- 
ing of the instabilities of the ground-state, i.e. nodeless, 
spatially localised solutions of Eqs. (1) under the action 
of the i-dependent perturbations. These solutions are 
well known, see, e.g. p7|-|3^ and refs. therein. Here, we 
restrict ourselves to the situation when the solitary waves 
are stable for dt = 0. Therefore we choose /3 > 0, because 
it ensures absence of the 'splitting' instability |30), and 
Vj_ = idx, to avoid collapse [[7pj34[|. 

It is important for the following to summarize relevant 
symmetry properties of Eqs. (1) with suppressed time 
derivatives {dt = 0). Invariance with respect to the two- 
parameter gauge transformation 

{E,,E2)^{Eie''^\E2e'*n, (2) 



leads to conservation of the energies Qi,2 = J dx\Ei,2\'^ 
or their equivalent combinations. There are also invari- 
ances with respect to transverse translation and Galilean 
transformation, 



£'1,2(2;) -^ £'1,2(2; + xo), 

£i,2(x) ^ £i,2(x - «z)e'''("-"^/2)^ 



(3) 
(4) 



01, 02, xq and V are free parameters. Although we use 
below the fact of the presence of the translational and 
Galilean symmetries, we do not need explicit expressions 
for the associated integrals of motion, which are linear 
momentum and 'center of mass', see, e.g., H. 

Symmetry property (g) indicates that the solitary so- 
lutions can be presented in the form 



£1,2(2;, z) = Ai^2{x)e' 



(5) 



^1,2(2;) are real functions obeying the system of ordinary 
differential equations 



5^^1,2 = Kl,2^1,2 



{Al^ + I3AIM 



1.2- 



(6) 



Exponential localization of the solitons requires ki,2 > 0. 
Actually one of these parameters can always be scaled 
away, which means that fixing one of them and varying 
the other in the whole region of the solitary wave exis- 
tence one will capture all possible situations. However, 
for convenience of analytical calculations it is better to 
keep them both. 

To study MI we seek solutions of Eq. (1) in the form 
of spatial solitons weakly modulated in time at frequency 

Ei,2{x,z) = {Ai,2{x) + 

{Ui.2{x, z) + iWi,2{x, z)) COS f7t)e'''i'^^+*'^ (7) 

Presenting solution of the linearized real system for the 
small perturbations t/m, Wm in the form Um ^ Um{x)e^^ 
and Wm ^ w,n{x)e^'^ we obtain the following eigenvalue 
problem (EVP) 



(£1 + -fn^i)u=^ -Aw, 

(£0 + -fn'^i)w = Xu, 



(8) 



where u — (ui, U2)"^, w — {wi, W2)'^ and / is the identity 
operator. Cq and £1 are: 



£0 



Vi 
Q 1)2 )' 



(9) 



£1 = 



-2(3AiA2 



-213A1A2 
B2 



(10) 



where I?i.: 



192 

2 X 



Ki,2 - ^1.2 - /3^2 1 and Si, 2 



By means of simple transformation one can reduce 
EVP (H) to the following two EVPs for real and imagi- 
nary parts of the perturbations: 



\2 ^ 

-A u. 



(£1 + in'^i){Co + 7f72/)w = -x^w. 



(11) 

(12) 

EVPs (O), (jl^) are adjoint to each other. Therefore 
they have identical spectra and in case of instability the 
imaginary and real parts of perturbations grow with the 
same rates. To answer stability question it is thus enough 
to study only one of the EVPs, and we concentrate below 
on the EVP (|l|). 

Let suppose that (Ki,2-f 7^^^) = ^1,2 > 0. Then, gener- 
ally, A^ € (— oOj ~^V) is a continuous part of the spectrum 
with unbounded eigcnfunctions, where Xg — min(^i,^2)- 
For particular cases when £1 becomes a diagonal oper- 
ator the continuum splits into two independent bands, 
(—00, —^i 2), corresponding to unboundness of ^1(2;) and 
^2(2;), respectively. Eigenvalues which do not belong to 
the continuum constitute the discrete part of the spec- 
trum and have bounded eigcnfunctions. Stable eigen- 
modes with eigenvalues obeying — Ag < A^ < 0, called 
'gap modes'. Any other mode of the discrete spectrum, 
i.e. any eigenmode with A^ complex or positive, renders 
the soliton unstable. If ^1 < and/or ^2 < 0, the gap is 
closed, Xg = Q. 

The procedure which we mainly follow to study stabil- 
ity of different types of solitary solutions consists of three 
basic steps. First, using analytical and numerical anal- 
ysis we identify the discrete spectrum of EVP ( p^ ) for 
ri = 0. Second, we develop perturbation theory for the 
neutral eigenmodes in the low- frequency limit, fi <C 1. 
Third, we numerically build continuations of all discrete 
eigenvalues into the region of finite positive Vi. We also 
allow for possible splitting of discrete eigenvalues from 
the edge of the continuum, but this was never actually 
observed. 



IV. INSTABILITIES OF CIRCULARLY 
POLARIZED AND MANAKOV SOLITONS 

The single wave solitons of Eqs. (ra) corresponding to 
right and left circular polarized e/m waves are 



Ai{x) — \/2Kisech\/2Kix, A2 = 0, 
Ai = 0, A2{x) = \j2K,2sech\J 2k,2X. 



(13) 
(14) 



For these solutions, EVP (gj|), separates into two inde- 
pendent scalar problems. Considering for example sta- 
bility of soliton with Ai ^ we get: 



(A^i + 7f72)(AAo + -i^^)w^ = -X^wu 

2 



-^dl + K2 + jn^-pAl 



W2 



\2 
-A W2, 



(15) 
(16) 



-^di + n,^2~iAi,~pAi,. 



where Hq = -\dl + ki - Al, Ni = -\dl + ki - 3Af . 



The operator on left-hand side of (|16|) has a non- 
negative spectrum therefore corresponding values of A^ 
(n = 0, 1, 2, 3 . . .) must be nonpositive, which means ab- 
sence of the unstable eigenmodes. In fact eigenvalue 
problem ( p^ can be solved analytically, see, e.g. p8| , p5[ . 
The eigenvalues are 



A^ = - f K2 + 7^' - Y [V8/3+1 - 2n - 1 



(17) 



Eq. (15) is exactly an EVP arising in theory of MI of 
solitons in single NLS equation and details of its analyti- 
cal and numerical investigations can be found in P,@,eI| . 
For the sake of completeness and comparison with MI 
of the other types of solutions we summarize the main 
results here. 

The discrete spectrum of operator A/iA/q consists of 
two neutral eigenmodes which can be readily identified 
by applying infinitesimal phase and Galilean transforms 
to solitary wave solution. These modes are wi^^ ~ A\ 
and w\v = xA\. Infinitesimal translations and variations 
of K\ generate two neutrally stable modes of the adjoint 
operator A/qA/i which are ui^ = dxAi and wi^^ = d/^-^Ai. 
These modes obey the identities: A/qWi^^ = 0, J\foWiy = 

-Uix, A/lUlKi_= -U'l0i, A/lMla- = 0. 

Following pit] we assume i7^ <C 1 and substitute the 
asymptotic expansions 



wi 



„(o) 



2,„(1) 



iw';" + n'w^'' + ...) 



and 



A2 = f72(A(i)' + r!2A(2)' + ...) 



(18) 



(19) 



into Eq. (fj). In the first two orders we have 



A/'iA/'ow;(°) = and ^fl^fow^^ 



^A(i 



^(0) 
(0) 

'j\ - 



li^fa 



A/i)u;(°^ Solution in leading order is w\"' — C^^wi^^ 
CvWiv, where C^^, Cy are constants. Orthogonality 
properties (wi0i,uia;) — (uiiujUiki) = (here and be- 
low (/,g) — J2m I d,^fmgm) rcsult in the independence 
of the branches produced by the phase and Galilean neu- 
tral modes. Therefore C^^ and C^ are in fact independent 
constants. Solvability condition of the 1st order problem 
gives 



i(i) 



^^^ =27. 



^(1)2 ^ r, {uix,Uix) 



47K1 



-^1-1- 



(20) 
(21) 



Eqs. (poD, ( |2l| ) indicate onset of instability for either sign 
of 7. However, the character of the instability depends on 
the sign of 7. For anomalous GVD (7 > 0) the spatially 
symmetric eigenmode becomes unstable leading to clus- 
tering of the soliton stripe into filaments (neck MI). For 
normal GVD (7 < 0) an excitation of the antisymmetric 
eigenmode leads to spatial symmetry breaking and bend- 
ing of the solitary stripe along the temporal coordinate 



(snake MI). The period of the modulations is approxi- 
mately equal to 2TT/njnaxj where ^max is the maximally 
unstable frequency. 

Typical dependencies of the MI growth rates vs fl are 
presented at Fig . H. The neck instability disappears at 
i^0i = ■\/3«;i/7 where wi = and ui = sech^^/2Kix. 
For 51 > ri^j the corresponding eigenmode becomes a 
gap mode. Note that for 7 > the gap becomes wider 
with increasing of fi. The snake instability disappear- 
ance is very difficult to track numerically because the 
corresponding eigenmode develops oscillating tails and 
becomes weakly localized, so that a larger number of the 
grid points is required. However, our numerical analysis 
clearly indicates that the branch of snake MI does not 
disappear stepwise at the point where the gap is closed, 
^3 = \/— K1/7, as was suggested in Ref. [pj, but contin- 
ues beyond this point and probably reaches A^ = at 
some larger fl. 

The nonlinear stage of MI is also perfectly analogous to 
that in the single NLS. Filaments, formed as the result of 
the neck MI development, collapse to a singularity during 
further propagation Rpl. The snake MI leads to soliton 
spreading due to unbalanced action of the self-focusing 
nonlinearity and normal GVD S . The second field E2 is 
not affected by the discussed instabilities, because of the 
incoherent nature of the coupling between Ei and £2- 

In the special case /3 = 1, ai = 0:2 and 71 — 72 Eqs.(l) 
are invariant under the arbitrary rotations in {Ei,E2) 
plane, Ei^2 -^ cost?i?i^2 ± sin^i?2,i. This leads to a new 
parameterization of the ground state solitons. These are 
usually called Manakov solitons Pq| and they are given 
by the solutions of Eqs. (pf) with ki_2 = i^- 



Ai = cosM(a;), A2 = sindA{x). 



(22) 



Here the angle is a new free parameter characterizing 
the polarization angle, and A{x) = \/2Ksech\/2Kx. Be- 
cause of the rotational invariance, Manakov solitons with 
different polarizations are equivalent and their MI is inde- 
pent of the polarization angle. Therefore one can always 
set 6 — 0, and then the corresponding EVP coincides 
with Eqs. (|l5l), ((lel). 



V. INSTABILITIES OF LINEARLY AND 
ELLIPTICALLY POLARIZED SOLITONS {(3 / 1) 

A. Soliton family and associated neutral modes 

To study solitons of an arbitrarily polarization for 
/3 ^ 1, i.e. ^1 7^ and A2 ^ 0, it is more convenient 
to introduce absolute, ip — ^{4>i + </>2), and relative, 
"0 = ^('^1 ~ 4'2), phases. The corresponding integrals 
of motion are the total energy Q = Qi + Q2 and energy 
unbalance Qu = Qi — Q2- Associated soliton parameters 
are k = ^(ki + K2) and S — ^{ki — K2). 

For S = there is an obvious and well known analytical 
solution of Eqs. (|) §,||-|| 



Ai^2ix) ^A{x) 



2k 



sechv2KX, 



(23) 



corresponding to a linearly polarized soliton. Using nu- 
merical solution of Eq. (o) one can verify that for /3 ^ 1 
the exact solution (|2^) belongs to the family of the soli- 
tary solutions parameterized by k and S. 

Using Eq. ( pj\) for n ^ 0, fl = and its analog for the 
solution (nj) we conclude that for fixed values of k and 
f3, a family of ground state coupled solitary solutions of 
Eqs. <M) exists for S G {—5c, 5c), where 



Sc 



l-4f5 + VT+W 



3 -f 4/3 - VT+W 



(24) 



Analogs of Eqs. ( p4| ) have been derived before in a num- 
ber of papers using different methods, see, e.g. [|9|, p7| , p2[ . 
Expression under the modulus in Eq. (24) changes its 



sign from plus to minus once (3 changes from /3 < 1 to 
/3 > 1. It follows that for /3 < 1 the family of the ellipti- 
cally polarized solitons splits from the family A2 = 
{Ai = 0) of the circularly polarized ones at 5 = 5c 
{5 = —5c) and this is vice versa for /3 > 1. Continuous 
variation of 5 from —5c to 6c for fixed k and f3 < 1 {j3 > 1) 
results in monotonic decay of Q2 (Qi) from its maxi- 
mal value Q+ (Q-) down to zero and in growth of Qi 
{Q2) from zero upto Q+ {Q~), where Q± — 2^2(k± 5c). 
Therefore, we can make an important for the following 
conclusion, that for f3 < 1 dsQu > and for /3 > 1 
dsQu < 0. Numerically build dependencies of Q„ vs 5 
for different values of (3 are presented in Fig. 0. 

Consider now the main spectral properties of the el- 
liptically polarized solitons for fi = and /3 7^ 0,1. 
Phase and Galilean symmetries generate three neutral 
eigenmodes of the EVP (|2|), they are w^ = (Ai,yl2)^, 
w^ = {Ai,—A2)'^, and w^ — x{Ai, A2)'^ . Infinitesimal 
variations of k and 5, and translational symmetry gen- 
erate neutral modes of the adjoint problem (pi]): u^ = 
d^{Ai,A2)^, us = 55(^1,^2)'^, and u^ = d^{Ai,A2)'^. 
These six modes obey the following identities CoWip — 0, 
Cqw^ = 0, CoWy = -Ux, CiUk = -w^, Cius = -w^, 
CiUx = 0. 

For /3 = Eqs. (1) separate in two independent NLS 
equations. The independence of the two fields results in 
additional translational and Galilean symmetries char- 
acterizing freedom of the relative transverse translation 
and motion of the two waves. Therefore EVPs ( [ll| ) and 
(O) have additional neutral modes usx — dx{Ai, —A2)'^, 
wsv = x{Ax, —A2)'^. As numerical solution for < /? < 
1 shows the corresponding eigenvalue produces stable 
branch of the discrete spectrum. For |/3| <C 1 approx- 
imate expression for this eigenvalue can be readily found 
[^,|o|, A^^ = — 64/J/15. Excitation of the corresponding 
eigenmode results in position oscillations of the soliton 
upon its propagation ||3l[]. When B — > 1 this eigenmode 
disappears into the continuum |3l|. 



B. Asymptotic stability analysis {Q.'^ <^ 1) 

Now assuming that /3 3> il^ we can use the asymptotic 
techniques described in the previous section to continue 
zero-eigenvalue modes into the region fi^ ^ 1. Making 
substitutions 



w = («?("' + n^w^-^'^ + ...) 



and of Eq. (|l|) 
ders: CiCqw^^^ 
7 (A - 



(25) 



into ([L2p we get in the first two or- 



and £i£ow(^) 



-Ad 



w 



(0) 



Ci)w^°\ Solution in the leading order is w*^°^ = 
CipW^p + Csws + CyWy. As in the previous subsection one 
can show that the branches produced by the two phase 
modes on the one hand and by the Galilean mode on the 
other are independent. The solvability condition of the 
first order problem for the Galilean mode gives 



A«' = -2r 






(26) 



which implies snake instability for 7 < 0. When 5 = 0, 



X\;' = -47K/3, cf. Eq. (|2 

For the two phase modes the solvability condition re- 
sults in a quadratic equation for A*-^^ 



where 



aA(i)' + 5A(i)' + c = 0, 



4a = d^QudsQ - d^QdsQr, 



(27) 



26 = jQid^Q + dsQu) - iQuid^Qu + dsQ), 



c = -Aj^QiQ2. 
Corresponding values of C^ and C^ are linked through 

flip pniifllifA;" 



the equality 



Cy _ 2-fQ-X(^')^dsQu 
C0 ~ X(^fdsQ-2-fQj 



where A'^' is the corresponding root of Eq. (|27| ). In the 
general case expressions for the roots of Eq. ( ]27| ) can not 
be analyzed analytically, but it is already clear that four 
symmetric neck type eigenmodes exist and either two or 
all of them may be responsible for instability. 



If ,5 = then 



,w^,us) 



{w^, 



and therefore 



w^ and Wtp eigenmodes produce independent branches 
of the discrete spectrum. This results in independence 
between C^ and C. 
dated eigenvalues: 



between C^ and C^ and simplifies formulas for the asso- 



A« =27 



d^Q 



dsQu 



47K, 



27k/(/3). 



(28) 
(29) 



Here /(/3) = (/ sechx g{x, (3) dx)~^ and function g obeys 
to {d1 — 1 + 2j^sech?x)g = sechx. f{j3) changes its 
sign from plus to minus when (3 passes through unity, 
see Fig. ^ Ahernatively, Eq. (^) can be rewritten as 



A 






2-fQi/dsQi, cf. Eq. (4) in || 



The A^ eigenvalue and associated neutral mode w^, 
are linked to the symmetry in the absolute phase ip and 
have their analogies in the spectral problem for single 
wave solitons described in the previous section, see Eq. 

(pO|). The AL^ eigenvalue and neutral mode w^ are 
novel. They can be directly attributed to the symmetry 
in the differential phase ip. This branch of the discrete 
spectrum generates instability for normal GVD (7 < 0) 
if /3 > 1 and for anomalous GVD if /3 < 1, see Figs. 
^, Q Thus, the asymptotic analysis indicates that for 
/3 > 1, 7 < neck and snake instabilities coexist, and for 
/3 < 1, 7 > two different types of neck instability coex- 
ist. Numerical evaluation of the roots of Eq. (E^ shows 
that the same conclusions hold also for S y^ 0, through- 
out the whole existence region of the family of cUiptically 
polarized solitons. Solving the EVP (|2|) numerically, we 
find that in the low-frequency limit the instability growth 
rates match those predicted by our perturbation theory 
within of few percent up to 17 ~ 0.5. Numerical investi- 
gation (for more details see below) also shows that apart 
from the three instabilities discussed in previous subsec- 
tion, a fourth MI associated with continuation of wsv into 
the region of /3 7^ 0, fi 7^ also exists. Analytical treat- 
ment of this instability is also possible, but it will not be 
pursued here, because corresponding MI branch is never 
dominant. 

Let us first discuss in general terms the physical mean- 
ing of all the different types of the instabilities in the 
simple situation with zero imbalancing [6 — 0), and only 
then we will proceed with details of the numerical anal- 
ysis. 



C. Instability induced polarization dynamics 



All eigenmodes of the EVP (12) are two component 
vectors, whose first and second components are respon- 
sible for the spatial form of modulations of the fields Ei 
and £'2, respectively. The eigenmodes w^p and w^ corre- 
sponding to the variations of the absolute phase ip and of 
the absolute velocity v of the coupled solitons have first 
and second components which are in phase for any value 
of X. This property holds throughout the whole region of 
existence of the associated branches of the discrete spec- 
trum. Therefore an excitation of these eigenmodes is not 
accompanied by the breaking of the polarization of the 
initial state. In contrast, the eigenmodes linked with w^ 
and wsv neutral modes, or in other words with varia- 
tions of the relative phase "0 and the relative velocity dv, 
have anti-phased first and second components. There- 
fore their excitation does lead to polarization symmetry 



breaking. In particular, one should expect that desta- 
bilization of the eigenmode associated with the relative 
phase ip will result in breaking of the linearly polarized 
soliton stripe into a chain of circularly polarized clusters, 
where neighboring clusters have opposite (left and right) 
polarizations. The same conclusions are obviously valid 
for the real parts of the perturbations. 

Note here, that by direct substitution of the linearly 
polarized solution Eq. (|2|) into the EVPs ([l|), (|l|) one 
can easily show that the eigenvalues of the eigenmodes 
with in-phase first and second components are indepen- 
dent of /3. This means physically that in-phase MI is in- 
sensitive to the relative strength of self- and cross-phase 
modulations. 



D. Numerical results for normal GVD (7 < 0) 

We start description of our numerical results with dis- 
cussion of the normal GVD case. We found two snake in- 
stabilities for /3 < 1, see Fig. 4. One of them corresponds 
to the in-phase snaking of both fields, see Fig. 4(c), and 
its growth rate in the low frequency limit is given by the 
Eq. (Pq). The other one corresponds to the anti-phase 
snaking, see Figs. 4(d). Examples of the growth rate 
dependencies vs fl and details of the anti-phase snaking 
appearance are presented in Fig. 4(a) and Fig. 5, respec- 
tively. Dependencies of the maximal instability growth 
rates on [3 are presented in Fig. H(b). We found that in- 
phase snaking dominates anti-phase snaking for all values 
of S and (3. For (3=1 the anti-phase snake mode dis- 
appears inside the continuum and does not appear again 
for aU /3 > 1. 

The dominant role of the in-phase snake instability 
means that breaking of the polarization state imposed by 
the initial conditions is unlikely to happen upon propaga- 
tion. Introducing imbalancing for a fixed total energy en- 
hances this dominance, see Fig. 4(a). Thus, when P < 1, 
the linearly polarized solitons are more stable than any 
other state of polarization. 

The anti-phase neck MI associated with the relative 
phase -0 appears for /3 > 1, see Eq. (29) and Figs. 6(b). 
The in-phase snake instability obviously also exists, see 
Eq. (|2^) and Figs. 6(c),(d). The in-phase snake MI 
dominates the anti-phase neck for 1 < f3 < j3sn and vice 
versa for [3 > (3sn, see Fig. 7, where the cross-over value 
f3sn depends weakly on S. This fact can also be seen 
from the comparison of the perturbative results for 6 = 0. 
According to the Eqs. (|8|), (||) and Fig. (0) the neck 
instability dominates the snake in the low frequency limit 
starting from P ~ 3.55. Numerical stability analysis gives 
that f3sn — 3.47 at fl = flmax for 6 = 0. Introducing 
imbalancing always leads to the suppression of the both 
instabilities, see Figs. 6(a), (c) and Fig. 7. Therefore 
the circular polarized soliton is most stable for a given 
energy. 

In analogy with MI of circularly polarized solitons for 



normal GVD, the neck and snake unstable eigenmodes 
become weakly confined and develop oscillating tails as 
n increases beyond the point where the gap is closed, 

To test our linear stability analysis and study the non- 
linear evolution we performed a series of computer simu- 
lations of the Eqs. (1) with initial conditions in the form 
of a soliton stripe perturbed by spatio-temporal white 
noise of order of few percent. Typical simulation results 
are presented in Fig. ||,y,|l^. For /3 < P^n we observed 
in-phase snaking of the stripe along the temporal dimen- 
sion, see Fig. ||. For (3 > (3sn the soliton stripe breaks 
in such a way as to form the interleaved intensity peaks 
of El and £'2, see Fig. g, as expected when the out of 
phase neck MI is dominant. The spatio-temporal pat- 
terns formed at the initial stage of MI finally spread be- 
cause of the unbalanced action of the normal GVD and 
self-focusing nonlinearity. For /3 ~ (3sn we observed com- 
petition between the neck and snake Mis, see Fig. |l^. In 
Fig. |l^ (bi), (b2) one can clearly see that at the inter- 
mediate stage of MI the typical in-phase snake pattern is 
superimposed on the anti-phase neck pattern. 

Thus, we conclude, that in the media with normal 
GVD spatial soliton stripes perfectly develop snake MI 
without polarization symmetry breaking if /3 e (0, I3sn) 
and neck MI with polarization symmetry breaking if 

/3e {f3sn,+00). 



E. Numerical results for anomalous GVD (7 > 0) 

There are two neck Mis in this case for /3 < 1, see 
Fig. nU One of them is associated with the absolute 
phase (fi and corresponds to the in-phase neck MI. The 
other one is associated with the relative phase tp and 
corresponds to the anti-phase neck MI. The in-phase MI 
dominates the anti-phase one for any value of S and (3, 
which means conservation of the polarization state im- 
posed by the initial conditions. Nonzero imbalancing for 
fixed total energy leads to the growth of the in-phase MI 
and to the suppression of the anti-phase one, see Fig. 
ll(a),(b). For f3 > 1 only in-phase instability exists, but 
now imbalancing leads to the suppression of the insta- 
bility, see Fig. 12. Presence of these instabilities agrees 
with the predictions of low-frequency analysis, see Eqs. 
(p8|), ( p^ ) and Figs. ^,||. A typical result of the numerical 
simulation of the neck instability development is shown 
in Fig. 1^. Note that the numerically attainable propa- 
gation distance was limited by the distance at which the 
most intense filaments formed at the initial stage of MI 
collapse to singularities. 

Cut-off frequencies, where the neck Mis disappear can 
be found analytically for (5 = 0. Growth rate of the in- 
phase MI becomes zero at 51 = -^8^/7 in full analogy 
with single NLS equation, see section IV. The anti-phase 
MI disappears at fi^ = k{D — B)/{2'-f) having w = and 
Ml 2 = (sec/iv^x)(^-i)/2^ here D = (11 - 5/3)/(l -I- (3), 



i3 = V(25-7/3)/(l + /3). 

Thus, in the media with anomalous GVD spatial soli- 
ton stripe always develops neck MI without polarization 
symmetry breaking and filamentary structure formed 
during this process collapses upon propagation. Detailed 
study of collapse in coupled NLS equations is outwith 
the scope of this paper. Some details on this issue can 
be found in IM. 



VI. DISCUSSION 

It is interesting to compare MI of solitons with re- 
sults on MI of CWs 1 13 1, which can be easily recov- 
ered from Eqs. (|), (]ll|), {^ putting d^ = 0. For 
simplicity we again consider the case of the linear po- 
larization, Ei,2 = \/k/(1 -t- /?)e*''^. Then correspond- 



ing eigenvalues are A^ = 757^ (2k — jil'^) and A?, 



i> 



7fl2(2K(l - /3)/(l + (3) - jfl^). For normal GVD, A^ 
can be positive only for /? > 1. For anomalous GVD, A^ 
generates instability for any f3 and A3, only for /3 < 1. 
Thus, as one could expect, neck instabilities of solitons 
related to the phase symmetries have analogies for CWs. 
Snake instabilities are obviously absent for CWs, which 
is the main difference between the dynamics of spatially 
confined solitons and infinitely extended CWs. Namely, 
in the case of normal GVD, CWs are modulationally sta- 
ble for /? < 1 and unstable for /3 > 1 (demonstrating 
polarization symmetry breaking). Solitons are snake un- 
stable in this situation for any /3 and this instability does 
not involve changes in the polarization state. However, 
starting from a critical value oi (3 — (3sn the snake in- 
stability becomes suppressed by the neck one, which is 
analogous to instability of CW. This instability does lead 
to polarization symmetry breaking. In particular, a lin- 
early polarized soliton breaks, due to this instability, into 
the chain of circularly polarized clusters. Because snake 
instability leads to spatial symmetry breaking and neck 
MI does not, the change in MI of solitons at /3 = (3sn 
can be interpreted as a transition from spatial symmetry 
breaking to polarization symmetry breaking. 

In the limit situation (3^1 self-phase effects are neg- 
ligible compare to cross-phase ones and development of 
the in-phasc and anti-phase neck Mis can be qualitatively 
explained using Fermat's principle. Due to MI develop- 
ment the effective refractive index for Ei and E2 fields 
gets modulated through the XPM mechanism with pe- 
riod 2Tr /rimax- This results in temporal cross-defocusing 
of filaments in media with normal GVD and in cross- 
focusing for anomalous GVD. Thus, in the case of nor- 
mal GVD, interleaved pattern of the intensity peaks of 
El and E2 fields should be preferable because it enables 
each field to see a refractive index that increases to its 
peak, i.e. one that is in accord with Fermat's principle. 
This is clearly verified in Fig. 9. The same arguments lead 
to the conclusion that a pattern with all intensity peaks 
coincident is preferable for anomalous GVD. 



Considering possibility of experimental observations of 
predicted phenomena, we have to say that diffraction in- 
duced MI of soliton-like stripe, which is formally equiva- 
lent to the case of anomalous GVD, is probably easiest to 
observe. However, it is less interesting at the same time 
because it is perfectly analogous to MI of CWs and it is 
not accompanied by any polarization effects. More inter- 
esting dynamics is expected in media with normal GVD. 
In fact, experimental observation of temporal splitting 
induced by normal GVD of spatially confined pulses in a 
self- focusing medium was recently reported in [ p7| . How- 
ever, transverse and polarization effects, which, accord- 
ingly to our results, should play an important role, were 
not studied during this experiment. Numerical studies 
IpTJ presented to support the experimental results were 
restricted by scalar approximation and radial geometry. 

The rescaled instability growth rate A as function of 
the modulational frequency fl can be related to physical 
units using the formulae: 



Aph — 



A 



4fci 



O^ - 

"p/i — 



7O2 



2kk"i 



Here Xph and Qph are the instability growth rate and 
modulational frequency in physical units, k is the wave 
vector, w is the beam width, k — d^k. For example 
for radiation at Ifim propagating in an AlGaAs planar 
waveguide k ~ —10~^^s^ / m ||38| and for typical soliton 
transverse size w ~ 50/i?n p9| we get Xph ~ A/ (5cm) 
and rjp^ ~ f7^/(10~^^s^). For experiments with fused 
silica at wavelength SSOrtm, see second from Refs. [p7[, 
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VII. SUMMARY 

We have analyzed and described dispersive MI of fam- 
ilies of nodeless spatial solitons in the system of the two 
incoherently coupled NLS equations. Considering cou- 
pled soliton states, we have established the existence of 
the four branches of instabilities, which are linked to the 
symmetries in the absolute and relative phases and in 
the absolute and relative motions of solitons. We gave 
a physical interpretation of our results describing GVD- 
induced polarization dynamics. In particularly, we found 
that in media with normal GVD the MI induced spa- 
tial symmetry breaking in the transverse plane changes 
to the polarization symmetry breaking when the relative 
strength of the cross-phase modulation exceeds a certain 
threshold value. In media with anomalous GVD, MI re- 
sults in breaking of spatial solitons into spatio-temporal 
clusters which collapse upon further propagation. This is 
not followed by either spatial or polarization symmetry 
breaking. 
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FIG. 2. Energy unbalancing Q„ vs 5, k = 1. 
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FIG. 3. Function /(/3), see Eq. (M). 



FIG. 1. Instability growth rates of the circularly polarized 
soliton vs O, Ki = 1. Full (dot-dashed) line is for neck (snake) 
MI, 7 = 0.5 (7 = -0.5). 
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FIG. 4. Instability growth rates, spatial profiles of the soli- 
tary solutions and of the unstable eigenmodes for /3 < 1, 
7 — —0.5. Dash-dot (dash-dot-dot-dot) lines correspond to 
the in-phase (anti-phase) snake MI. (a) Growth rates vs fl, 
13 = 0.3. Thin (thick) lines correspond to n = 1, 5 = Q, 
Q ~ 4.04, Q„ = (k = 1.155, 5 = 0.5, Q ~ 4.04, Q„ ~ 2.18). 
(b) Maximal growth rate vs /3. Thin (thick) lines correspond 
to 5 — {S = 0.5). (c) Components of the eigenmode 
corresponding to the in-phase snake MI, /3 — 0.3, S = 0, 
n = 1. (d) Components of the eigenmode corresponding to 
the anti-phase snake MI, f3 = 0.3, 5 = 0, n = 0.92. 
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FIG. 6. Instability growth rates, spatial profiles of the soli- 
tary solutions and of the unstable eigenfunctions for /3 > 1, 
7 = —0.5. Dash-dot (dash) lines correspond to in-phase snake 
(anti-phase neck) Mis. (a) Growth rates vs fi, /3 = 2. Thin 
(thick) lines correspond to k = 1, 5 = 0, Q = 1.75, Qv. = 
(k = 0.93, <5 = 0.2, Q = 1.75, Q„ = -0.91). (b) Compo- 
nents of the eigenmode corresponding to in-phase snake MI, 
/3 = 2, 5 = 0.2, f2 = 0.8. (c) Growth rates vs fi, /3 = 7. Thin 
(thick) lines correspond to k = 1, 5 = 0, Q ~ 0.68, Qv. = 
(k = 0.97, (5 = 0.2, Q ~ 0.68, Q„ ~ -0.16). (d) Components 
of the eigenmode corresponding to anti-phase neck MI, /3 = 7, 
(5 = 0.2, n = 1. 
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FIG. 5. Eigenvalues corresponding to the anti-phase snake 
MI vs fl for several choices of /3: n = 1, 5 = 0, -y = —0.5. 
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FIG. 7. Maximal growth rates of the in-phase snake 
(dash-dot line) and the anti-phase neck (dash line) Mis vs 
j3 for /? > 1, 7 = —0.5. Thin (thick) lines correspond to 5 = 
{& = 0.5). Cross-over occurs at Pan ^ 3.47 for (5 = and 
Pan ~ 4.04 for 5 = 0.5. 
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FIG. 8. Development of the in-phase snake MI for 13 = 2, 
K = 1, 5 = 0, 7 = -0.5. (ai,2) |£i,2| for z = 12; (bi,2) |Si.2 
for z = 14.7. 
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FIG. 9. Development of the anti-phase neck MI for /3 = 7, 
K = 1, 5 = 0, 7 = -0.5. (ai,2) IS1.2I for z = 8.4; (bi,2) l^i.al 
for z = 10.2; (01,2) |£i,2| for 2 = 12.6. 
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FIG. 10. Competition between the in-phase snake and 
anti-phase neck Mis: (5 = 3.47, k = 1, 5 = 0, 7 = —0.5. 
(ai,2) \Eia\ for 2 = 9; (bi,2) |Si,2| for z = 12 



FIG. 11. Instability growth rates, spatial profiles of the 
solitary solutions and of the unstable eigenfunctions for /3 < 1, 
7 = 0.5. Full (dash) lines correspond to in-phase (anti-phase) 
neck Mis. (a) Growth rates vs f2, /3 = 0.3. Thin (thick) lines 
correspond to k = 1, (5 = 0, Q ~ 4.04, Qu = (k = 1.155, 
5 = 0.5, Q ~ 4.04, Qu ~ 2.18). (b) Maximal growth rate vs 
13. Thin (thick) lines correspond Xo 5 = G {5 = 0.5). (c) Com- 
ponents of the eigenmode corresponding to in-phase snake 
MI, P = 0.3, 5 = G, Q. = 1.5. (d) Componets of the eigen- 
mode corresponding to anti-phase snake MI, /3 = 0.3, 5 = 0, 
= 0.9. 
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FIG. 12. In-phase MI growth rates for /3 > 1, 7 = 0.5. 
(a) Growth rates vs O, /3 = 2. Thin (thick) lines correspond 
to K = 1, 5 = 0, ~ 4.04, (5„ = (k = 1.155, 5 = 0.5, 
Q ~ 4.04, Qu ~ 2.18). (b) Maximal growth rate vs /3. Thin 
(thick) hues correspond to 5 = (5 = 0.5). 
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FIG. 13. Development of the in-phase neck MI for /3 = 2, 
K = 1, 5 = 0, 7 = 0.5. (ai,2) 1^1,2 1 for z = 4.6. With further 
increasing of z most intense filaments develop collapse. 
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